#delimit;
clear all;
version 13.1;
pause on;
program drop _all;
capture log close;
set more off;

/*INSERT FOLDER PATH*/
cd /;

/**/


foreach defn in strct {;

use multi_level_dataset_005_field_vars.dta, clear;
drop if ovrlp_anypos==1;
keep id newid pmid year nbcites_it rltd rtrct_year score treat age srce_pmid pubyear case_code yr subv*;

compress;
gen yr_to_rtrct=year-rtrct_year;

forvalues i=5(-1)1 {;
	gen yrs_before_rtrct`i'=0;
	replace yrs_before_rtrct`i'=1 if yr_to_rtrct==-`i' & rltd==1;
};


forvalues i=1(1)10 {;
	gen yrs_after_rtrct`i'=0;
	replace yrs_after_rtrct`i'=1 if yr_to_rtrct==`i' & rltd==1;
};
gen yrs_after_rtrct11=0;
replace yrs_after_rtrct11=1 if yr_to_rtrct>10 & rltd==1;

gen yrs_before_rtrct6=0;
replace yrs_before_rtrct6=1 if yr_to_rtrct<-5 & rltd==1;

drop yr_to_rtrct;
order id newid pmid srce_pmid rltd pubyear year nbcites_it rtrct_year score treat age yr case_code subv_of_trth_strct subv_of_trth_loose yrs_before_rtrct6;

xtqmlp nbcites_it yrs_before_rtrct6-yrs_before_rtrct1 yrs_after_rtrct1-yrs_after_rtrct11 i.age i.yr if subv_of_trth_`defn'==1, fe i(id) cluster(case_code);
estimates save estimates/poisson_main_dynamics_`defn'.ster, replace;

estimates drop _all;
estimates use estimates/poisson_main_dynamics_`defn'.ster;

matrix B=e(b);
matrix V=e(V);

keep if rltd==1;
gen t_to_rtrct=year-rtrct_year;
bysort t_to_rtrct: keep if _n==_N;

foreach y in 1 2 3 4 5 {;
	local t=colnumb(B,"yrs_before_rtrct`y'");
	gen yrs_before_rtrct`y'_coef = B[1,`t'];
	gen yrs_before_rtrct`y'_var = V[`t',`t'];
};

foreach y in 1 2 3 4 5 6 7 8 9 10 {;
	local t=colnumb(B,"yrs_after_rtrct`y'");
	gen yrs_after_rtrct`y'_coef = B[1,`t'];
	gen yrs_after_rtrct`y'_var = V[`t',`t'];
};

gen treatment_coefs=0;
gen treatment_vars=0;

foreach var in 	yrs_before_rtrct5 yrs_before_rtrct4 yrs_before_rtrct3 yrs_before_rtrct2 yrs_before_rtrct1 
 		yrs_after_rtrct1 yrs_after_rtrct2 yrs_after_rtrct3 yrs_after_rtrct4 yrs_after_rtrct5 
 		yrs_after_rtrct6 yrs_after_rtrct7 yrs_after_rtrct8 yrs_after_rtrct9 yrs_after_rtrct10 
 		{;
		replace treatment_coefs=`var'_coef  if `var'==1;
		replace treatment_vars=`var'_var  if `var'==1;
		};

gen treatment_95low=treatment_coefs-1.96*sqrt(treatment_vars);
gen treatment_95hi=treatment_coefs+1.96*sqrt(treatment_vars);

sort t_to_rtrct;

twoway scatteri -0.25 0 0.25 0, c(l) lcolor(black) lwidth(thin) msym(none) yline(0, lcolor(black) lwidth(thin))
|| line treatment_95hi t_to_rtrct, lpattern(shortdash) lcolor(red) lwidth(medium) || line treatment_95low t_to_rtrct, 
lpattern(shortdash) lcolor(red) lwidth(medium) || line treatment_coefs t_to_rtrct, lpattern(solid) lcolor(blue) lwidth(medthick) || 
if t_to_rtrct >= -5 & t_to_rtrct <= 10, xtitle(" " "Time to Retraction") xlabel(-5(1)10, labsize(small)) ylabel(-0.25(0.05)0.25, 
format(%15.2f) angle(horizontal) labsize(small)) legend(off) saving(graphs/poisson_main_dynamics_`defn'.gph, replace);

twoway scatteri -0.25 0 0.25 0, c(l) lcolor(black) lwidth(thin) msym(none) yline(0, lcolor(black) lwidth(thin)) || rcap treatment_95hi treatment_95low t_to_rtrct, lcolor(blue) || scatter treatment_coefs t_to_rtrct, msize(medsmall) mcolor(red) || if t_to_rtrct >= -5 & t_to_rtrct <= 10, xtitle(" " "Time to Retraction") xlabel(-5(1)10, labsize(small)) ylabel(-0.25(0.05)0.25, format(%15.2f) angle(horizontal) labsize(small)) legend(off) saving(graphs/poisson_main_dynamics_`defn'_rcap.gph, replace);

};

/*************/


/*for (figure 4) rank interactions*/

foreach defn in loose strct {;

use multi_level_dataset_005_field_vars.dta, clear;
drop if ovrlp_anypos==1;
keep if subv_of_trth_`defn'==1;
keep id newid pmid year nbcites_it rltd rtrct_year ranking treat age srce_pmid pubyear srce_pmid yr case_code subv*;
compress;


forvalues i = 1(5)100 {;
                display `i';
                local j=`i'+4;
                gen byte treatXrnk`i' = treat*(ranking==`i'|ranking==`i'+1|ranking==`i'+2|ranking==`i'+3|ranking==`i'+4);
                label variable treatXrnk`i' "After Retraction � Rltd'ness Rank btw. `i' & `j'";
                };

gen byte treatXrnk100plus = treat*(ranking>100);
xtqmlp nbcites_it treatXrnk1-treatXrnk96 treatXrnk100plus i.age i.yr, fe i(id) cluster(case_code);
estimates save estimates/poisson_main_rnk_intr_`defn'.ster, replace;

estimates drop _all;
estimates use estimates/poisson_main_rnk_intr_`defn'.ster;

matrix B=e(b);
matrix V=e(V);

bysort ranking: keep if _n==_N;

foreach y of numlist 1 6 11 16 21 26 31 36 41 46 51 56 61 66 71 76 81 86 91 96 {;
	local t=colnumb(B,"treatXrnk`y'");
	gen qntl_score`y'_coef = B[1,`t'];
	gen qntl_score`y'_var = V[`t',`t'];
};

gen qntl_score100plus_coef=B[1,colnumb(B,"treatXrnk100plus")];
gen qntl_score100plus_var=V[colnumb(B,"treatXrnk100plus"),colnumb(B,"treatXrnk100plus")];

gen qntl_coefs=.;
gen qntl_vars=.;

foreach j of numlist 1 6 11 16 21 26 31 36 41 46 51 56 61 66 71 76 81 86 91 96 {;
	replace qntl_coefs=qntl_score`j'_coef  if ranking==`j';
	replace qntl_vars =qntl_score`j'_var   if ranking==`j';
};

replace qntl_coefs=qntl_score100plus_coef if ranking>100;
replace qntl_vars=qntl_score100plus_var if ranking>100;

gen qntl_95lo=qntl_coefs-1.96*sqrt(qntl_vars);
gen qntl_95hi=qntl_coefs+1.96*sqrt(qntl_vars);

sort ranking;

keep if inlist(ranking, 1,6,11,16,21,26,31,36,41,46,51,56,61,66,71,76,81,86,91,96,101);
keep ranking qntl_95lo qntl_95hi qntl_coefs;
replace ranking=ranking-1 if ranking>1;

gen rkng=string(ranking[_n])+"-"+string(ranking[_n+1]);
replace rkng="100+" if ranking==100;
sencode rkng, replace;

twoway rcap qntl_95hi qntl_95lo rkng || scatter qntl_coefs rkng, msize(medsmall) mcolor(teal) yline(0, lcolor(black) lwidth(thin)) || if ranking >= 1 & ranking <= 100, ytitle(" " "Magnitude of the Treatment Effect", size(medsmall)) xtitle(" " "Relatedness Rank") ylabel(-0.50(0.10)0.50, format(%15.2f) angle(horizontal) labsize(small)) xlabel(#21, valuelabels labsize(vsmall) labgap(*2) angle(45) tlength(*.5)) legend(off) note(" " "Rank is {it:decreasing} in relatedness", span size(vsmall)) saving(graphs/poisson_main_rnk_intr_`defn'.gph, replace);

};